Self heating and nonlinear current- voltage characteristics in bilayer graphene 
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We demonstrate by experiments and numerical simulations that the low-temperature current- 
voltage characteristics in diffusive bilayer graphene (BLG) exhibit a strong superlinearity at finite 
bias voltages. The superlinearity is weakly dependent on doping and on the length of the graphene 
sample. This effect can be understood as a result of Joule heating. It is stronger in BLG than in 
monolayer graphene (MLG), since the conductivity of BLG is more sensitive to temperature due to 
the higher density of electronic states at the Dirac point. 
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I. INTRODUCTION 

Two-dimensional graphene in its monolayer and bi- 
layer forms can exhibit rather different electronic 
characteristics.SS In monolayer graphene (MLG) the 
valence and conductance bands touch each other at 
two inequivalent Dirac points in the Brillouin zone, 
around which the bands are linear. Thus the density of 
states (DOS) also vanishes linearly around these points, 
where the Fermi energy of charge-neutral graphene is lo- 
cated. In the most common (Bernal) stacking of bilayer 
graphene (BLG) the two layers are electronically coupled 
such that the two linear bands of the individual layers 
mix to form four bands, the lower of which are parabolic 
around the Dirac points. In this case the DOS around 
these points is approximately constant. 

This difference gives rise to different charge screening 
and transport properties for the two types of graphene. 
In particular, the temperature dependences for the con- 
ductivity of diffusive MLG and BLG around the charge- 
neutrality point (CNP) differ.lSHll In both cases, thermal 
excitation of quasiparticles from the valence to the con- 
duction band (i.e. thermal creation of electron-like and 
hole-like charge carriers) increases the conductivity with 
temperature, as is typical for semiconductors and insu- 
lators. However, due to differences in the DOS, in MLG 
the conductivity at CNP grows only quadratically with 
temperature, while in BLG it increases linearly.^S 

In this paper we show how this difference of MLG and 
BLG is reflected in the current-voltage {I{V)) character- 
istics of diffusive graphene. For MLG it is known that 
the I{V) characteristics tend to be linear at low bias volt- 
age V and have a tendency to saturate at higher volt ages 
due to scattering of electrons from optical phonons.l^ESl 
Close to CNP the I{V) at small V can become super- 
linear as a result of Zener-Klein tunneling between the 
valence and conduction bands, especially in low-mobility 
samples.'Sl By measuring the I{V) curves of both MLG 
and BLG on a Si02 substrate in a two-terminal con- 
figuration, we show that in BLG the I{V) characteris- 
tics have a much stronger tendency for superlinearity at 
V < 0.1 Y, which we associate with an increase of the 



conductivity due to self heating (Joule heating) . This ef- 
fect is only weakly dependent on the level of doping and 
on the length of the sample. We confirm this interpre- 
tation with numerical simulations using a semiclassical 
model based on Boltzmann theory for a diffusive two- 
dimensional (2D) system in quasiequilibrium. The model 
takes into account electron scattering from charged im- 
purities, the band-bending effects due to charge doping 
by the metallic source and drain electrodes,'^^"^ uniform 
impurity doping, as well as nonuniform doping by a gate 
electrode. For MLG the Joule-heating-related nonlinear- 
ity is found to be weak, which is consistent with previous 
experiments and calculations .ISHn] 

The paper is organized as follows. In Sec. [ll] we in- 
troduce the theoretical model. Sec. |III| describes the ex- 
perimental results and compares them to theory, and in 
Sec.|TV]we end with some discussion of other mechanisms 
for current nonlinearities. Details of the model are given 
in the Appendixes. In App. |A] solution of the electro- 
static part of the problem in terms of a Green function 
is detailed. Appendix |B] discusses the modeling of the 
impurity scattering and gives expressions for the charge 
density and the transport coefficients for MLG and BLG. 
Finally, in Sec. [C] we discuss analytic semiclassical results 
for the temperature dependence of the conductance of a 
p-n junction, which can form in the graphene close to the 
metallic electrodes. 

II. THEORETICAL MODEL 

The model geometry that we consider (see Fig. [ij is a 
simplification of typical experimental geometries. It con- 
sists of a two-dimensional (quasi-one-dimensional) chan- 
nel coupled to two transport electrodes {I and r) as well as 
a back gate electrode (bg) . In this geometry we solve for 
the electrostatic potential together with transport equa- 
tions for the quasiparticles in the channel. 

Our semiclassical transport theory assumes a diffusive 
2D electron system in quasiequilibrium j'i^ such that the 
quasiparticle distribution is described by a local chemical 
potential fi{x) and a local temperature T{x). A third 
unknown field is the electrostatic potential f{x) in the 
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FIG. 1. The rectangular geometry considered in the model. I 
and r depict the left and right (source and drain) electrodes, 
while bg is a back gate electrode. All three are held at different 
constant potentials, (j>h 4'r, and (fjbg, respectively. The 2D 
channel, of length L, is at a distance d from the back gate. 
For simplicity the regions below and above the channel are 
assumed to be occupied by the same dielectric medium with 
permittivity e. 



channel {y — d), which shifts the energy Ed{x) of the 
local Dirac point such that Eo{x) — —e(p{x). The three 
fields (p{x), n{x), and T{x) are solved from the equations 

nL/2 

ip{x) = / d^G{x, d; ^, d)e[n{^) - ndop]/e 

J-L/2 

+ X! '^x{x,d:)(f>x (1) 

X=Lr,bg 

[aix)^i'{x)/e + -f{x)rix)y = 0, 
-[a{x)n\x)/e + k{x)T'{x)]' = Pj{x). 

Here prime denotes the x-derivative, G is the Green func- 
tion of the Laplace operator, ipxix, y) is the "characteris- 
tic function" for electrode X (see App.[A|, and e — £r£o, 
where £r and eq are the relative and vacuum permittivi- 
ties. Further, n{x) is the two-dimensional charge density 
in units of the electron charge — e (App.[B|), and Udop is a 
phenomenological doping density that describes the dop- 
ing effect by impuriti es, wh ich is assumed to be constant. 
(Thus, charge puddle^^^EZl are not described — see below 
for discussion). The factors a{x) and k{x) are the charge 
and thermal conductivities, respectively, while a{x) and 
7(x) are the thermoelectric transport coefficients, satisfy- 
ing a{x) = T{x)^{x). These factors are inversely propor- 
tional to the impurity density riimp, which is taken to be 
another constant parameter independent of Udop (App. 
[B| . The quantities n{x)^ cr(x), k{x), ^{x), and a{x) all de- 
pend on i^{x). Finally, Pj{x) = j{x){fi' {x)/e) is the Joule 
power per area, with j{x) = u{x)^' {x) / e + ^(x)T'{x) 
being the electric current density, which is conserved, 
j{x) = j. The total current through a system of finite 
width > L, d is / = jW. 

The boundary conditions at x = ±L/2 are chosen to 
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FIG. 2. Results for BLG (a,b) and MLG (c,d) with parameter 
L = 350 nm, fi^g = 0.05 eV, Udop = 10.0 • 10^'' m'^ Sr = 4.0, 
d = 210 nm. Additionally for BLG n^^^p = 7.0- 10^*^ m"^, and 
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(a,c) Linear-response conduc- 



tivity a(V = 0) = {L/W){dl /dV)v^o at three temperatures 
and (b,d) the corresponding finite-V^ differential conductiv- 
ity a{V) = {L/W){dI/dV) for To = 4 K at indicated gate 
voltages (j)bg. 



which correspond to a voltage bias V, assuming negli- 
gible contact resistance, the right-hand electrode to re- 
main grounded, and both electrodes to remain at the 
bath temperature Tq. The finite equilibrium chemical 
potential fieq 7^ describes the effects of the work func- 
tion mismatch and the resulting charge transfer between 
the graphene and the electrodes, with fi^^ > {ii^q < 0) 
leading to n-type (p-type) doping.^** 

The first of Eq. ([T]) is the Poisson equation written as 
an integral equation at y = d. The second and third 
describe current conservation and heat balance, respec- 
tively. Note that in order to keep the model simple, we do 
not include electron-phonon coupling and thus the Joule 
power is dissipated only via diffusion. This should be 
reasonable first approximation at bias voltages V well 
below optical phonon energies.G^ 

In Fig. [2] we show typical results obtained from the 
model for BLG and MLG with parameters obtained from 
the fit to (the gate dependence of) our BLG experiments 
in Sec. |III[ It is to be noted that the aspect ratio L/d = 
1.7 is far too small for a simplified parallel-plate capacitor 
model to work properly in the geometry of Fig. [l] and a 
full numerical solution of Eqs. ([T]) is needed. The fields 
obtained as their solutions are very nonuniform, as shown 
for BLG in Fig. [3| 

The upper panels of Fig. [2] are for the case of BLG. Fig. 
[2|[a) shows the gate dependence of the linear-response 

conductivity a{V = 0) = {L~^ j p{x')dx')^^ , where 

p(x) — l/a{x) is the resistivity. Due to finite positive 
doping density {ridop > 0) the point of minimal conduc- 
tivity (the apparent "CNP") is shifted to = 4>bg,min ~ 
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FIG. 3. Bias dependence of the field profiles for n-type BLG 
at gate voltage (j)tg = 8 V: (a) conductivity a, (b) temperature 
T (c) charge density n (d) chemical potential fi and the local 
Dirac point En (higher and lower curve, respectively). Solid, 
dashed, and dash-dotted curves are for V = —0.1, 0.0, and 
0.1 V, respectively. Dotted lines show additionally a, n, and 
Ed in equilibrium at (j>hg — 0; here n{x = 0) « ridop- The 
parameters are as in Fig. |2] and we have defined the units 
o"o = CblgEo and no — (e/d)Eo/e^ , where Eq =1 eV (for 
Cblg, see App.|B|. 



-13 V and as result of the lead-doping effect (/igg > 0) 



the gate dependence is asymmetrical around the CNP.li^l 
Since the lead-doping is of n type, at (ptg > 4'bg,min the 
BLG is of n type everywhere. However, for < 4>hg,min 
two p-n junctions appear j^i^l and the BLG becomes of n-p- 
n type. While for (pjjg > 4>bg,min the dependence on tem- 
perature is relatively weak, for (p^g w 4>bg,min it is roughly 
linear: (7{V = 0) oc T (Ref. For (j)bg < (f)bg,min, 
where the p-n junctions dominate the resistance, it can 
be shown that approximately a{V = 0) cx — l/ln(r) as 
T — >■ (App.[C|. However, the theory is valid only when 
all parts of graphene remain far from the local CNP. This 
is because charge puddles, quantum-mechanical effects 
(see App. [C]), and a possible gap in the BLG spectrum 
are not taken into account. Thus, we mainly concentrate 
on the gate voltages 0bg > 4>bg,m,in- 

Fig. I2|b) shows the differential conductivity (t{V) — 
{L/W){dI/dV) = L{dj/dV) as a function of V for a 
few gate voltages (j)bg > 4>bg,min at low bath tempera- 
ture. Note that the (j{V) curves are nearly symmetric 
[a{~V) « <y{V)]^ although there are small deviations, 
which are due to the asymmetrical choice of the bound- 
ary conditions and the presence of the gate electrode. 
The increase of (y{V) at finite V signifies a superlinear 
contribution to the I{V) curve: I{V) « dV + 6*2^^ 
iy > 0) with G2 > 0. This superhnearity is strongest 
close to the CNP, where the conductivity of BLG is most 
sensitive to temperature (see Eq. ( |B5[ )). In fact, the in- 
crease of the conductivity with voltage [see Fig. |3][a)] is 
entirely due to heating, which leads to maximal temper- 
atures of T{x) ~ 300 K at y = 0.1 V [Fig. ^h)]. The 



charge density n{x) for example, remains almost indepen- 
dent of V [Fig. [sj^c)]. Indeed, the density is of the form 
n{x) = h{fi{x) — Ed{x)) (App.jBj), and n{x) — E£i{x) re- 
mains close to its value at ^ = everywhere [Fig. [3|^d)]. 
Thus the bias voltage "gates" the graphene very little. 
The dotted lines in Fig. |3] additionally show the results 
at equilibrium, with (pbg = 0. The fast transients close to 
the electrodes are due to the doping by the leads. This 
doping is not restricted only to the region on the order of 
a screening length 1 nm (Sec. |b| from the leads, but 
is actually long-ranged.'SEIl 



The lower panels Fig. [2][c,d) show equivalent results 
for MLG, where the impurity density nimp has been cho- 
sen so that the conductivities are of a similar magnitude 
as for BLG. The temperature-dependence of a{V = 0) 
for (pbg ^ 4>bg,min is uow clcarly even weaker. For <j)bg ~ 
4'bg,min it is quadratic, a{V — 0) cc (Ref. [51), and for 
't'bg < 4>bg,min hucar, a{V = 0) cx T (App. [C|. Corre- 
spondingly, the increase of the cr{V) at (f>bg >q)bg,min is 
much weaker. This is consistent with the fact that the 
I{V) curves measured for MLG are typically linear or 
even sublinear, except close to CNP in low-mobility sam- 
ples, where Zener-Klein tunneling is of importance.'^ In 
the case of MLG the "gating" effect of the bias voltage is 
somewhat larger due to the longer screening length, but 
remains also weak. 



Here we have concentrated on short samples, with 
L/d ^ 1. In the considered model geometry the gate 
dopes the graphene quite weakly at distances on the 
order of d from the ends. Additionally, the center of 
the graphene heats more than the ends. Therefore at 
4>bg > 4>bg,min the euds tend to dominate the resistance 
[Fig. [3]ja)]. When L ^ d, the parallel-plate limit is ap- 
proached, where n{x) and (t{x) become uniform, with 
fj-{x) and Ei:i{x) roughly linear. One may then sim- 
plify the equations Eqs. ([I]) by taking a = 7 = 
and using the Wiedemann-Franz law k — CTa, where 
C — (7r^/3)(fcB/e)^, and assuming a constant a. The 
temperature profile is thus approximated with T{x) — 
y/T§ + [1/4 - {x/L)^]V^/C, which scales simply with V. 
The heating effect on the conductivity (j{V) therefore de- 
pends relatively weakly on the length of the sample. 



We do not pursue further simplifications or extensions 
of the model here, but it should be noted that the ther- 
moelectric coefficients a and 7 are not of great impor- 
tance for the current nonlinearity. However, under some 
conditions the strong temperature gradient at the ends 
can also cause the conductivity to decrease at small bias 
voltage. A very weak sign of this is seen in the flat region 
of the (j)bg = 18 V curve in Fig. [2]jb). We also note that 
our tests with some simple models for charge puddles can 
reduce the width of this flat region, making nonlinearity 
stronger also at high gate voltages. 
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FIG. 4. (a) Optical picture of two bilayer graphene samples 
of lengths L = 350 nm and L = 950 nm. The BLG flake 
has been colored in red to enhance its visibility, (b) Raman 
spectrum of the corresponding bilayer graphene sheet 
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FIG. 5. Typical I{V) curves for MLG and BLG at the CNP 
at 4.2 K. The dashed line has a slope equal to the conductivity 
4.5eV^ of both samples at V = 0. For MLG the /(V) curve 
is linear. In BLG it is superlinear, which we associate with 
the Joule heating of the sample. 



III. EXPERIMENTS 

We have measured seven bilayer graphene samples in a 
two-lead configuration and found qualitatively the same 
transport properties for each of them. Here we focus 
on the results obtained on two samples from the same 
BLG sheet having lengths L = 350 nm and L = 950 nm, 
and widths W = 900 nm and 1550 nm, respectively (Fig. 
|4|. The samples were contacted using Ti/Al/Ti sandwich 
structures with thicknesses 10 nm / 70 nm / 5 nm (10 nm 
of Ti is the contact layer). Three 0.6 /im wide contacts 
were patterned using e-beam lithography. The strongly 
doped Si substrate was used as a back-gate, separated by 
270 nm of Si02 from the sample. 

The I{V) curve of our 0.95 x 1.55 ^m^ sample at 4.2 
K is illustrated in Fig. [5] together with the I{V) in a typ- 
ical MLG sample. While the MLG result is linear,I^Il2l 
the BLG curve exhibits clearly superlinear behavior at 
small drain-source voltage V; the nonlinearity in BLG 
depends relatively weakly on the gate voltage (fibg but 
is largest near the CNP. The nonlinear behavior can be 
observed more clearly in Fig. iG^a) which displays the dif- 
ferential conductivity a{V) =~(L/W){dI /dV) of the long 
and short samples at a few values of (l)bg- It is seen that 
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FIG. 6. The left-hand panel (a) shows the experimentally 
measured differential conductivity cr{V) at finite V for both 
the short (dashed line) and the long (solid line) BLG sample at 
three pairs of gate voltages that are chosen so that the minima 
at F = V for both samples coincide. The right-hand panel 
(b) shows the corresponding zero-bias conductivity a{V = 0) 
vs. the gate voltage (/)6g for both samples. The temperature 
is To = 4.2 K. 



also the length dependence of the nonlinearity at bias 
voltages below « 0.1 V is weak. This supports its in- 
terpretation as a heating effect: as mentioned above, the 
temperature should scale with V and not, for example, 
the electric field V/L. Note, furthermore, that the cr{V) 
curves are very symmetrical. 

Figure [6]jb) shows the full gate voltage dependence of 
the zero-bias conductivity of the short and long sam- 
ples. In both samples, the minimal conductivity is lo- 
cated in the negative gate voltage region around -10 
V. The minimum zero-bias conductivities are roughly 
3.8 and 4.7e^//i for the short and long sample, respec- 
tively. These are close to the value ^ Ae^ /h typically 
found for both MLG and BLG.'^'^^' An asymmetry be- 
tween the n-doped and p-doped regions is clearly visi- 
ble and is more pronounced for the short sample where 
the conductivity is almost constant in the p-region. We 
interpret this electron-hole a symm etry as a sign of the 
leads n doping the graphene j'^^'^ so that there are p-n 
junctions present at larger negative gate voltages. This 
is consistent with expectations for Ti/Al electrodes.'^^nil 
In a parallel-plate approximation the charge density is 

~ {Cbg/e){(j)bg-(l)bg,min), where Cbg = 1.3-10"'' F/m^. 
Using this we estimate from the slope of aiV = 0) vs. 
n > for the long sample the mobility fj,„i — cr/(en) to 
be at least 1500 cm^V"^s~^. Using this the mean free 



path is estimated to be l„ifp — ^imh^m/ e < 30 nm, and 
thus the samples are diffusive. 

We compare the raw experimental data to the theoret- 
ical results in Fig. [Tj where solid and dashed lines are for 
experimental and theoretical data, respectively. Figures 
^a) and (b) show the a{V = 0) vs. 4>bg dependence for 
the short and the long sample, respectively. Also exper- 
imental data measured at 77 K and 300 K are shown. It 
is seen that the zero-bias conductivity increases slightly 
from 4 K to 77 K, and more significantly from 77 K to 300 
K. The conductivity change is strongest near the CNP 
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FIG. 7. Experimental results (solid lines) and a theoretical 
fit (dashed lines) for a short BLG sample (L = 350 nm) and a 
long BLG sample {L = 950 nm). Left-hand panels (a,b) show 
the linear response conductivity at T = 300 K (top curve), 
77 K, and 4 K (bottom curve), and the right-hand panel (c) 
shows the differential conductivity for L = 350 nm at To = 4 
K and at gate voltages 8 V (top curve), 3 V, and —5 V (bottom 
curve). The parameters used for the theoretical fit are: Heq = 
0.05 eV, riir^p = 7.0-10^® m"^, Sr = 4.0, d = 210 nm. For the 

while 

2 



'Pbg ^ (f'bg.min, where some parts of the system are close 
to the CNP. At these gate voltages the low-temperature 
theoretical results for (t{V — 0) tend to fall well below 
the experimental ones, whereas at 300 K the opposite is 
true. In addition to this the clearest discrepancy is that 
the self-heating predicted by the theory is too strong, as 
shown by the overly steep slope of (t{V) in Fig.[7][c). Also, 
the slight decrease of the experimental slope with voltage 
is not captured. A proper explanation of these effects 
would require considering interactions of the electrons 
with phonons, partic ularly the remote interface phonons 
of the Si02 substrate.G^^EH Another clear deviation is that 
the theoretical (y{V) curves for high 0{,g tend to be flatter 
close to V — than in the experiments. (The theoretical 
results for the long sample are otherwise similar as in Fig. 
7[c), but even shghtly more "flat".) As suggested in Sec. 
n] the correct shape could presumably be reproduced by 
considering charge puddles, i.e., a non- uniform impurity 
density and doping. To keep the model simple and the 
number of fitting parameters small, we have neglected 
puddles here. 

It should be noted the gross features of the results 
can be understood simply based on the temperature- 
dependence of the local conductivity,^ which may be 



short sample L = 350 nm and Udop = 10.0 • 10 ^ m 
for the long sample L = 950 nm, and ridop ~ 7.0 • 10^^ m 



worked out analytically [Eq. (B5) below] and the fact 
that the average temperature increases roughly linearly 
with the bias voltage. However, the precise shape of 
the nonlinearity is dependent on the various sources of 
nonuniformity. 



and also in the p-doped region. Since the theory is ex- 
pected to work only in the absence of p-n junctions, the 
fitting is done only to the a{V — 0) vs. (ptg data [Fig. 
[7]ja,b)] at (pbg > 4'bg,min, with /igg > 0. The same param- 
eters are then used for calculating cr{V). The a{V) data 
for the short sample are shown in Fig.[7]jc). 

The parameters used for the fit are given in the cap- 
tion of Fig. [Tj The work function mismatch fieq — 0.05 
eV is of the correct sign and order expected for Ti/Al 
electrodes .fSl A gate distance d = 210 nm is used, which 
is smaller than the experimental 270 nm. The smaller d 
used for the comparison is reasonable, since in the sim- 
plified geometry assumed by the theoretical model the 
gate potential is more strongly screened by the trans- 
port electrodes: even for the long sample the parallel 
plate limit is not fully reached.'^ It is also notable that 
the average doping density n^op = 7-10 • 10^^ 1/m-^ is 
much smaller than the impurity density riimp = 7.0 • 10^^ 
1/m^, although the simplest theories for impurity doping 
predict these two to be equal.l^ Even though we model, 
for simplicity, the impurities to be all similar and of the 
screened- Coulomb type (App.pL in reality there may be 
several different types of disordei=^ present in our samples, 
which complicates the situation. The orders of magni- 
tude ndop,nimp ~ 10^^-10^^ 1/m^, are still in the range 
that can be expected from other experiments.!^ 

The overall agreement of the theory with the ex- 
periment is good, apart from the large deviations for 



IV. DISCUSSION 

The heating effects described above are not the only 
possible sources of nonlinearity. As already mentioned, 
at high enough bias voltage electron scattering from 
phonons tends to reduce the conductivity, which is ex- 
pected to make the current-voltage curves at high bias 
sublinear. This has been seen in M LG as a tendency for 
the current to almost saturate.'^ESl We also see similar ef- 
fects in our experiments^ with BLG at voltages V > 0.1 
V. 

The possibility of nonlinear I{V)s in graphene at low 
bias have also been discussed based on simple argu- 
ments involving the energy-dependence of the number 
of open transport channels in the Landauer-Biittiker de- 
scription of the linear-response conductivity.^"^ Such 
calculations are problematic for the prediction of current- 
voltage characteristics beyond linear response, however, 
since they do not consider the role of the actual voltage 
profiles.'^ Ideally, the electrostatic potential drop should 
be calculated self consistently, as we have done above. 

Other mechanisms for nonlinear (mostly superlinear) 
current-voltage responses in grap hene h ave very recently 
been discussed by many authors.'^^^^lHlI] particular, 
Zener-Klein tunneling in MLG has been sh own to give 
rise to superlinearities close to CNP.EHMZ! It seems un- 
likely, however, that a similar mechanism would be of 
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importance in our experiments, since the nonlinearity 
is weakly gate-dependent. Other possibihties for non- 
hnearities include the presence of tunnel junction^^H or 
contact phenomenal^ but we can disregard them as an 
explanation of our measurements due to the weak length- 
dependence of the nonlinear conductivity cr{V). Fur- 
thermore, nonlinear current-voltage curves in graphene 
oxides have been explained with space-charge limited 
currents,!^ but such effects are likely to be negligible 
in metallic graphene, as also supported by the absence 
of any bias-doping effects in our simulations. Nonlincari- 
ties are predicted also for vertical transport in misaligned 
BLG or few-layer systems.!^ However, of all these pos- 
sibilities the self-heating scenario presented above seems 
to be the most likely explanation of our experimental re- 
sults. 

To summarize, we have shown how Joule heating can 
contribute to the shape of the observed current-voltage 
characteristics of bilayer graphene in the diffusive limit. 
The heating is signified by a strong superlinear contribu- 
tion in I{V) and thus a low-bias differential conductivity 
(t{V) increasing with V. Our experimental results and 
our numerical calculations are in good overall agreement 
for bias voltages < 0.1 V. The heating effect is much 



stronger in bilayer graphene than in monolayer, as can 
be expected from the differences in their electronic struc- 
tures. 
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Appendix A: Green function and characteristic 
functions 

The Green function G{x,y;x' ,y') of the Laplace op- 
erator satisfies \7^G{x,y; x' ,y') — 5{x — x')5{y — y'), 
G{x' ,y' ;x,y) — G{x,y;x' ,y') and zero boundary con- 
ditions on the electrodes. For the particular geometry 
under consideration (Fig. [ij the Green function may be 
found analytically using conformal mapping techniques 
and it is 



^ , 1 J / (sinicoshy — sini' coshy')^ -f (cosisinhy — cosi' sinhy')^ 

' ' ' 27r y (sin i cosh y — sin i' cosh y')^ -f (cos i sinhy -f- cos ai' sinhy')^. 



(Al) 



where x — ■kx/ L, y — ny/L and so on. 

The characteristic function i/jx for electrode X is de- 
fined such that it satisfies the Laplace equation V^Vx = 
with the boundary condition that Tpx = 1 on electrode 
X and ipx = on the other electrodes. Once the Green 
function is known, the characteristic functions may be 
easily represented in terms of it: 



f°° d 

ij,{x,y)= dy' — G{x,y;x',y')\^,^^^^ 
i^bgix,y)^- I dx' — G{x,y;x\y')\y,^g. 



(A2) 



-L/2 



Appendix B: Charge concentration and transport 
coefficients for graphene 

In our semiclassical model the charge density is as- 
sumed to be of the form 



(Bl) 



n{x) = / dEV^{E, x) [f{E - fi{x), T{x)) 
- 6{Ed{x) - E)], 



where f{E,T) = l/(e^/'=«^ -f 1) is the Fermi function, 
Ed{x) = —e(p{x), and we define Vip{E,x) = 'D{E + 
e(p{x)), where !?(£') is the density of states (DOS). Us- 
ing Eq. (Bl) self consistently in the Poisson equation is 
essentially the Thomas-Fermi approximation.'-^ 

Assuming a diffusive system with only elastic impurity 
scattering the transport coefficients are given by 



aix)= J dEa^{E,x)F{E - fi{x),T{x)) 
n{x)J-^ f dEa^{E,x)^p4S^PiE-^,{x),T{x)) 



7(x) 



ksTix) 



^ / dE<j^{E,x f J!^''^ F{E-fiix),Tix)), 
e J kbT^x) 



(B2) 



with a{x) = T{x)-f{x). Here F{E,T) ^ -df{E,T)/dE 
is the thermal broadening function, and we define 
a^iE,x) = a{E + e(p{x)), where a{E) = e^D{E)V{E) 
is the energy-dependent Boltzmann conductivity. The 
quantity D{E) — v'^{E)t{E)/2 the diffusion constant, 
where v{E) is the group velocity and t{E) the transport 
relaxation time. 

It is easy to see by change of the integration variable 
that the quantities in Eqs. (Bl I and (B2 1 only depend on 
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the difference fi{x) — Ed{x). Thus below we define the 
"hatted" quantities with n{x) = h{p{x) — Eo{x),T(x)), 
a{x) — (j{pL{x) — E£i{x),T{x)), and similarly for the other 
transport coefficients. 

As a specific model for the impurity scattering we 
consider only screened Coulomb impurities, which lead 
to a linear dependence of the conductivity on charge 
densitj0 for both BLG and MLG, as observed in most 
experiments!^ (For BLG also short-range scattering may 
be of importanceP) For simplicity we assume all of the 
bare impurities to carry a charge ±e and to be at zero 
distance from the graphene, and perform an average with 
respect to their positions.*^ For BLG and MLG some fur- 
ther approximations are made, as explained below. 



1. Bilayer graphene 

For BLG we assume a purely parabolic and gapless dis- 
persion E = ±{hvok)'^ /ji, where vq = 10^ m/s and 71 = 
0.4 eV. This yields a group velocity v{E) = 2voy/\E\/y[ 
and a constant DOS V{E) = hjrf^- Then 



_ 1 71 

TT (hVo)' 



(B3) 



The DOS leads to an inverse Thomas-Fermi screening 
length qTF,BLG = 2e^li/[4:Tre{hvo)^] ~ 1 nm"^ 

For the charged impurity scattering we use the 
"complete-screening" approximation.^^ Thus we find 
t{E) = ^ JJmyT^^ where riimp is the average impurity 
density. This yields 



a{E)^CBLG\E\, Cblg 
Then the transport coefficients are 
u{^i,T) =CBLG^kBTlTi (2 cosh 



8e^7i 



(B4) 



2kBT 



k(^,T) ^CTCBLG^kBThip/ksT) 



(B5) 



7(M,T) =CBLG'ikBT[- 

+ Li2(-e 
^2 1,2 



kBT 

-p./kBT\ ^ 



In (2 cosh( 



2kBT 



)) 



Li2(- 



where C — is the Lorenz number. Here and Li2 is 

the dilogarithm function. '^^ and h{a) = h{—a) is defined 
as 



h{a) = / \x\{x--af 



d 



1 



dx " + 1 



dx. 



(B6) 



This function has the limits h{a) « ^2 ln(2 cosh(a/2)), 
when |a| ::g> 1, and h{a) « 9^(3), when \a\ <C 1. Using 
these we see that the Wiedemann-Franz law k = CTa 
only applies if \^\ ^ fc^T. 



2. Monolayer graphene 

For MLG the dispersion relation is, E ^ zthvok, giving 
a constant group velocity v{E) — vq = 10^ m/s, and a 
density of states V{E) = ^ (Hfp • Then 



1 2{kBTf 
n{fi,T) = - ^^^^^^ giji/kBT). 



(B7) 



Here we have defined the function g{a) = Li2(— e^") — 
Li2(— e") — —g{—a), which has the limits g{a) « a|a|/2, 
when \a\ ^ 1, and g{a) ~ 21n(2)a, when \a\ <C 
1. The inverse screening length is now qtf.mlg — 
ikpe^ /{4:TTehvo), with kp = \fj,\/{hvo). 

For the impurity scattering we now assume that the 
"effective fine-structure constant"'^- of MLG is small, 
rs = QTFMLG/i'ikF) = e^/i^nehvo) < 1. (For Si02 
e = 4.0, and 0.5.) In this way we find t{E) — 

— - These give 



nimp 7r2 (hvo)^ r2 



'jiE) - C. 



MLG 



E\ 



C. 



MLG 



(B8) 



The transport coefficients are thus 



a{^i,T) =C 



MLG 



-{ksTY 



k{^i,T) ^CTCmlg 
2^2 



(kBT)' 



(B9) 



-CMLGfJ-kBT. 



Note again that the Wiedemann-Franz law k = CTa is 
only approximately valid in the limit 3> ksT. 



Appendix C: Low-bias resistance of p-n junction: 
classical thermal activation vs. quantum tunneling 

In order to understand the temperature-dependence of 
the conductivity in Fig. [2] at (j)f,g <^ (f>bg,min, we discuss 
some analytical results for the semiclassical conductance 
of a p-n junction in BLG or MLG. The existence of a p-n 
junction at location x — xq means that /leq — Eu^xq) = 
0. At low temperature we linearize E]j(x) around this 
point, such that ^eq — Eoix) ~ —Ai(x — xq), where Ai — 
E'jj{xq). The classical linear-response conductance for 
width W is then G — W[j_^p{x)dx)^^ , where p{x) — 



[a{Heq - Ed{x),T)] 



BLG 



In this case a{pL,T) is given by Eq. (B5|. Since we 



use the parabolic-band approximation, the x integral di- 
verges logarithmically and a cutoff length Lc is needed, 
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which should be on the order of L. In this way, the con- 
ductance of a p-n junction (width W) may be approxi- 
mated with 



G 



2\n{AiLj2kBTy 



(CI) 



The temperature-dependence has a logarithmic singular- 
ity at T = 0. This is the behavior seen in Fig. [2ja) at 

ft^bg ^ 4'bg,min- 

Clearly the semiclassical result must break down at 
low enough temperature, in which case some quantum- 
mechanical result taking into account Zener-Klein tun- 
neling is needed. The zero-temperature conductance 
would then remain finite. The simplest way to approx- 
imate the crossover temperature is to use a Wenzel- 
Kramers-Brillouin (WKB) approximation in a similar 
fashion as done for MLG.^^ '^^ Estimates of this type 
show that the crossover temperature may well be on the 
order of room temperature. We note that such a calcula- 



tion predicts a superlinear current / cx V"" with a = 4/3, 
unlike in MLG where a = 3/2. 



MLG 



Here a{ix, T) is found from Eq. (B9). The p-n junction 
in MLG has a conductance 



G = 



2AkBTCMLGAiW 



(C2) 



The linear temperature dependence is seen in Fig. [2][c) 
at (pbg <C 4>bg.min- At low temperature the p — n junc- 
tions completely dominate the conductivity of the en- 
tire sample. However, again, at low enough temperature 
this result breaks down. WKB estimates shows that this 
may occur already close to room temperature. Thus the 
Boltzmann calculations are only valid in the absence of 
p-n junctions. 
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